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A NUMERICAL  VISCOUS- INVISCID  INTERACTION 
ALGORITHM  FOR  PREDICTION  OF  NEAR-FIELD 
VSTOL  JET-INDUCED  FLOWFIELUS 
By 

A.  J.  Baker  and  P.  D.  Manhardt 
Computational  Mechanics  Consultants,  Inc. 
3601A  Chapman  Highway 
Knoxville,  TN  37920 


SUMMARY 


A viscous-inviscid  interaction  algorithm  has  been  derived  for  numerical 
determination  of  the  detailed  three-dimensional  flowfield  surrounding  and 
entrained  into  the  wake  of  a V/STOL  jet  issued  orthogonally  into  a subsonic 
crossflow  bounded  on  one  plane  by  an  aerodynamic  surface.  The  dominant 
influence  parameter  coupling  the  three-dimensional  external  potential  flow 
analysis,  and  the  imbedded  three-dimensional  viscous  and  turbulent  analysis  of 
the  jet  and  nearfield-wake,  is  the  entrainment  flow  induced  into  the  jet  region 
and  its  influence  on  the  potential  flow.  The  two  distinct  solution  domains 
are  pressure  coupled,  and  the  entrainment-induced  mass  flux  addition  to  the 
wake  is  the  primary  iteration  variable  of  the  global  algorithm. 

A key  element  of  the  algorithm  is  definition  of  a surface  defining  the 
coincident  boundary  of  the  dual  analyses.  Numerical  studies  have  confirmed 
existence  of  the  interface  surface  for  the  test  case  of  a circular  jet  issued 
into  a crossflow  with  velocity  ratio  equal  to  eight.  This  determination 
yields  definition  of  the  specific  geometry  of  the  interface  intersection  with 
the  injection  plane  and  the  entrainment  velocity  distribution  thereon.  Avail- 
able experimental  surface  pressure  data  are  specifically  employed  in  this 
analysis  phase,  and  available  additional  data  are  used  to  complete  definition 
of  the  interface  surface  evolution.  The  summary  three-dimensional  jet  region 
computations  yield  determination  of  the  detailed  evolution  of  the  turbulent 
flow  and  employs  a turbulence  closure  model  involving  all  six  components  of  the 
Reynolds  stress  tensor.  Finite  element  numerics  have  been  employed  to  compu- 
tationally evaluate  the  five  distinct  computational  phases  of  the  algorithm. 
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INTRODUCTION 


In  the  transition  from  hover  to  wing-borne  flight,  a significant 
fraction  of  the  lift  of  a V/STOL  aircraft  is  furnished  by  direct  engine 
thrust.  This  injection  of  a high-velocity  jet  of  exhaust  (or  air)  at  almost 
right  angles  to  the  aerodynamic  surface  and  to  the  cross-flow  induced  by 
the  aircraft  forward  flight,  produces  an  extremely  complicated  three- 
dimensional  local  flowfieid  which  can  significantly  alter  the  aircraft  aero- 
dynamic characteristics.  Since  the  basic  phenomena  of  transverse  jet 
injection  is  so  fundamental  to  V/STOL  performance,  it  is  crucially  important 
that  a firm  understanding  of  cause  and  effect  be  established.  The  mission 
of  this  research  project  is  to  derive,  implement  and  computationally  evaluate 
the  key  features  of  a three-dimensional  viscous-inviscid  flow  interaction 
numerical  solution  algorithm  for  the  near-field  of  a V/STOL  jet. 

Since  the  transverse  injected  jet  is  fundamental  to  many  physical 
processes,  in  addition  to  V/STOL  aircraft,  a comprehensive  base  of  experi- 
mental data  and  linearized  theoretical  analyses  has  become  established. 

Perhaps  the  first  consequential  theoretical  study  was  reported  in  the  disser- 
tation of  Chang  (ref.  1),  who  employed  potential  flow  theory  and  bound  vortex 
filament  concepts  to  predict  roll-up  of  the  seoaration  boundary  between  a 
parallel  onset  flow  and  a perpendicular  cylindrical  jet.  A series  expansion 
solution  yielded  the  prediction  of  the  initial  roll-up  near  the  jet  exit 
plane  illustrated  in  Figure  1.  A (hand)  numerical  solution  predicted  the 
further-field  roll-up  shown  in  Figure  2,  which  clearly  illustrates  evolution 
to  the  hallmark  horse-shoe  cross-sectional  shape.  This  basic  concept 
(application  of  potential  flow  theory)  has  been  expanded  upon  and  refined  to 
a great  extent,  and  has  produced  (viscous-corrected)  models  of  the  V/STOL  jet 
that  can  accurately  predict  gross  jet  flow  parameters  such  as  jet  centerline 
trajectory,  lateral  spread,  cross-section  shape  and  mean  axial  velocity  at 
each  cross-section  (c.f.,  ref.  2,  3). 

The  earliest  analyses  showed  only  the  potential  flow  aspect  of  the 
interaction  of  the  two  flows  while  current  models  employ  empirical  corrections 
to  approximate  the  dominant  viscous  effects  such  as  "entrainment."  Jordinson 
(ref.  4)  conducted  early  experiments  on  a round  jet  that  confirmed  the  horse- 
shoe cross-section  contours,  and  also  noted  that  the  jet  drew  the  flow  boundary 
layer  up  into  the  jet  wake  region.  This  is  not.  a dominant  feature  of  a jet 
issued  from  an  isolated  pipe  (ref.  5),  and  has  become  recognized  as  perhaps 
a most  consequential  contributor  to  performance  alterations  of  V/STOL  aircraft. 
Entrainment  is  properly  interpreted  as  the  total  effect  of  vortex  roll-up  plus 
the  induced  turbulent  mixing,  while  blockage  is  simply  the  characteristics  for 
inviscid  flow  around  an  equivalent  solid  body.  The  potential  flow  model 
equivalent  is  therefore  a cylinder  with  suction  coupled  with  an  empirical 
accounting  for  viscous  effects. 

The  significant  interest  in  the  V/STOL  jet  problem  prompted  the  conducting 
of  many  comprehensive  experimental  test  projects  to  probe  the  associated  three- 
dimensional  velocity  fields  to  well  downstream  of  injection  and  to  record 
induced  surface  pressure  distributions.  It  is  this  later  phenomenon  that  is 
perhaps  the  singularly  important  analysis  requirement  since  in  practice  the 
injection  plate  becomes  a (lifting)  aerodynamic  surface.  Vogler  (ref.  6), 
Bradbury  and  Wood,  (ref.  7)  published  early  results.  McMahan  and  Mosher 
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(ref.  8)  published  oil  flow  streaklines,  see  Figure  3,  that  provide  dramatic 
evidence  of  the  induced  flowfield  difference  between  a V/STOL  jet  and  the 
equivalent-diameter  solid  blockage.  For  the  latter.  Figure  3a,  definite 
stagnation  points  exist  both  upstream  and  downstream  on  the  center-line  and 
are  connected  by  an  apparent  inviscid  flow  streamline.  Outside  this  stream- 
line, the  exterior  flow  apparently  remains  irrotational  and  simply  responds 
to  the  blockage  induced  by  the  body.  The  flow  interior  to  this  streamline 
curves  inwards  to  the  centerline  and  appears  to  form  a viscous-dominated  closed 
recirculation  region.  Since  the  plate  is  impervious,  a (small)  axial  flow  must 
become  induced,  directed  away  from  the  plate  and  parallel  to  the  solid  cylinder 
for  some  distance.  Thereafter,  the  alternate  shedding  of  the  Karman  vortex 
street  probably  results  (c.f.,  ref.  5). 

The  plate  streamlines  for  the  circular  V/STOL  jet.  Figure  3b,  are  con- 
sequentially different  (note:  the  concentric  circle  is  part  of  the  experimental 
apparatus  and  not  a streakline).  The  appearance  of  the  upstream  stagnation 
point  region  does  not  appear  consequentially  different  indicating  blockage- 
dominance.  However,  in  the  mid-chord  region  of  the  jet,  the  surface  flow 
becomes  directed  almost  radially  inwards,  in  distinct  contrast  to  solid 
blockage  flow.  Two  streamlines  are  symmetrically  oriented  downstream,  which 
separate  the  complex  wake  interior  flow  from  the  deflected  irrotational  free- 
stream.  The  incoming  streamlines,  which  divide  the  flow  region  into  entrained 
or  deflected  segments,  intersect  the  downstream-dividing  streamlines  symmetrically 
and  a weak  stagnation  line  appears  to  connect  these  two  points  (essentially  on 
top  of  the  circle).  The  close-up  view  of  a similar  test  is  shown  in  Figure  4 
(from  ref.  9).  A rear  stagnation  point,  which  must  be  induced  by  flow  impinging 
the  plate  from  above  (anti-parallel  to  the  jet  axis)  is  clearly  evident,  as  is 
the  trajectory  of  a streamline  therefrom  directed  to  the  dividing  streamline  of 
the  onset  flow.  Interior  to  this  envelope,  the  flow  appears  directed  radially 
inwards  towards  the  jet,  which  from  conservation  of  mass  must  turn  parallel  to 
the  jet  and  form  the  characteristic  wake.  It  is  quite  apparent,  therefore, 
that  blockage  dominance  for  the  V/STOL  jet  is  confined  to  a small  upstream 
region,  and  that  the  majority  of  the  near-jet  surface  flow  is  highly  rotational, 
turbulent  and  three-dimensional.  Therefore,  one  might  expect  that  potential 
flow  models  would  encounter  considerable  difficulty  in  addressing  the  near  flow 
of  a V/STOL  jet. 

An  additional  consequential  experimental  data  base  has  been  generated  by 
probing  the  three-dimensional  flowfield  downstream  of  transverse  injection  for 
a circular  jet.  Mosher  (ref.  10)  concludes  that,  for  a circular  jet  with 
velocity  ratio  4 <_  x <_  10,  which  for  incompressible  flow  is  basically  the  ratio 
of  the  jet  mean  velocity  and  the  onset  reference  freestream  velocity,  entrain- 
ment is  the  primary  bending  mechanism  and  becomes  more  dominant  as  x is 
increased.  Pressure  bending  is  apparently  important  only  in  the  first  few  jet 
diameters  downstream  from  the  injection  plane.  He  further  established  that 
the  jet  path,  defined  by  the  location  of  the  counter-rotating  vortes  flow,  was 
a dominant  influence  on  the  near-jet  surface  pressure  distribution.  Referring 
to  Figure  5,  which  is  an  interpretation  of  the  jet  roll-up  process,  and  estab- 
lishment of  the  counter-rotating  vortex  pair,  Mosher  probed  extensively  the 
jet  flowfield  vicinity  in  three-dimensions  for  the  circular  jet  of  velocity 
ratio  x = 8.  However,  equipment  limitations  prevented  the  taking 


of  data  closer  to  the  plate  than  about  six  jet  diameters.  Kamotani  and 
Greber  (ref.  11  ) also  probed  the  field  for  A = 8 for  a circular  jet.  They 
pursued  considerable  data  interpretation,  in  particular  determination  of  the 
plane  orientation  along  the  jet  path  with  unit  normal  parallel  to  the  local 
jet  dominant  axial  velocity.  Plotting  data  in  vertical  planes  produced  the 
familiar  horse-shoe  profiles  previously  determined  by  Jordinson.  However, 
when  the  local  velocity  vector  is  resolved  into  scalar  components  in  the 
local  normal-transverse  plane  coordinate  system,  the  local  jet  cross-section 
shape  and  transverse  plane  counter-rotating  flow  immediately  becomes  apparent. 
Figure  6 illustrates  the  symmetric  half-plane  distribution  of  the  normal  mean 
axial  jet  flow  for  A = 8,  at  downstream  stations  xx/d  = 7 and  23,  where  xj 
is  the  curvilinear  coordinate  parallel  to  the  curved  jet  path  with  origin  at 
the  plate  surface,  and  d is  the  jet  initial  diameter.  Note  the  characteristic 
horse-shoe  shape  has  become  the  (equally)  familiar  "kidney"  shape  in  this 
representation.  Since  the  jet  is  basically  axially  symmetric  at  the  point  of 
injection,  the  combined  influence  of  blockage  plus  entrainment  is  to  rapidly 
erode  the  jet  potential  core  and  to  eventually  establish  the  axial  velocity 
maxima  near  the  center  of  each  transverse  plane  vortex,  Figure  6b. 

Figure  7 illustrates  the  projection  of  the  measured  velocity  field  at 
xj/d  = 23  onto  the  transverse  plane.  The  rotational  structure  is  clearly 
evident  as  is  existence  of  an  upstream  stagnation  point  in  the  plane. 
Insufficient  downstream  data  were  taken  to  resolve  transverse  plane  closure, 
but  an  inviscid  flow  separation  streamline  is  evidenced  in  this  plane.  The 
experiments  of  Fearn  and  Weston  (ref . 12  ) did  include  a sufficiently  large 
region,  and  Figure  8 illustrates  the  measured  transverse  plane  resolution  of 
the  velocity  vector  with  axial  velocity  isovels  superimposed  for  A = 8,  at 
axial  stations  x^d  = 5.2  and  8.3.  Both  upstream  and  downstream  stagnation 
points  are  evidenced  in  transverse  velocity,  and  they  appear  connected  by  an 
inviscid  flow  streamline  in  this  plane.  Therefore,  while  the  oil  flow 
streaklines  differ  appreciable  between  the  transverse  jet  and  solid  blockage 
experiments  for  A = 8,  it  appears  that  at  least  by  x:/d  = 5,  the  jet  surface 

has  become  inviscidly  impervious.  Alternatively  stated,  the  entrainment  of 
fluid  into  the  jet  wake  appears  primarily  completed  within  the  first  few  jet 
diameters,  and  thereafter  as  the  jet  bends  and  proceeds  downstream,  the  action 
of  viscosity  and  turbulence  homogenizes  the  flow.  Thus,  entrainment  appears 
primarily  the  result  of  the  existence  of  the  plate.  That  the  region  of  flow 
surrounding  the  jet  eventually  achieves  a downstream  transverse  plane  velocity 
stagnation  point  is  also  confirmed  by  the  oil  flow  streaklines  in  Figures  3-4. 
This  is  a plausible  mechanism  for  driving  the  flow  anti -parallel  to  the  jet 
(towards  the  plate)  behind  the  immediate  wake. 

Fearn  and  Weston  (ref.  13)  document  detailed  results  of  an  experimental 
project  to  determine  surface  pressure  distributions  induced  by  a 0.1m  (4-inch) 
diameter  circular  jet  over  the  velocity  ratio  range  2 £ A < 10.  Also  included 
are  specific  comparisons  to  O.N.E.R.A.  and  Mosher  data,  and"  the  results  are 
presented  as  detailed  Cp  distributions  on  radial  rays.  Their  results  show 
that  a region  of  positive  pressure  always  exists  on  the  upstream  side  of  jet. 
Large  radial  pressure  gradients  are  created  laterally  by  the  induced  mid-side 
low  pressures,  and  the  minimum  pressure  rotates  upstream  as  A increases. 
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Figure  7.  Transverse  Plane  Velocity  ufc  Vector  Plot  (Ref.  11)  For 
Circular  Jet,  A = 8,  Xj/d  = 23 
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The  radial  pressure  gradient  along  the  rearward  ray  is  typically  much  smaller, 
and  x appears  the  controlling  parameter. 

The  conclusions  and  indications  for  the  round-jet  experiments  appear 
to  carry  over  to  non-circular  jets  in  essence.  Mosher  (ref . 10)  presents 
oil  flow  streaklines  for  an  approximate  rectangular  cross-section  (with 
rounded  corners)  with  long  dimensions  both  parallel  and  perpendicular  to 
the  onset  flow.  As  expected,  the  cross-stream  orientation  produces  a 
relatively  larger  region  of  disturbed  flow  and  low  surface  pressures. 

Weston  and  Thames  (ref.  14  ) measured  the  plate  pressure  distributions  shown 
in  Figure  9 for  a rectangular  jet  of  aspect  ratio  five.  It  is  evident  that 
all  cases  exhibit  large  regions  of  induced  negative  pressures  with  extent  a 
function  of  aspect  ratio  and  orientation.  However,  it  seems  fair  to  conclude 
that  the  essential  interaction  between  blockage  and  entrainment  is  unchanged 
by  modest  jet  initial  geometry  modifications  for  the  velocity  ratios  evaluated. 

These  experimental  results  seem  to  indicate  that  the  non-potential 
character  of  the  near  jet  flow  is  the  dominant  determinor  of  the  essential 
V/STOL  jet  character.  Therefore,  it  seems  rather  optimistic  to  assume  that 
a viscous  correction  to  a potential  flow-derived  model  could  adequately 
describe  near-field  phenomena  without  an  intractible  number  of  adjustable 
parameters.  This  is  the  base  assumption  of  this  analysis  as  has  been  con- 
cluded by  others.  Specifically,  numerical  solutions  of  some  form  of  viscous, 
rotational  flow  governing  partial  differential  equations  is  required. 

Kotansky  and  Bower  (ref.  15)  present  a concise  summary  of  the  basic  dilema 
in  this  area  which  is  repeated  here  as  Table  1.  A quick  scanning  of  this 
table,  published  in  1975,  shows  that  the  numerical  technology  required  to 
predict  entrainment  and  cross-winds  is  long-term,  see  columns  5 and  6. 

The  alternative  is  correlation  of  results  from  a large  data  base  which  can 
be  quite  cumbersome  and  expensive,  see  Columns  1-2.  Kotansky  et.  al  (ref. 16  ) 
presents  a solution  approach  for  a full  V/STOL  aircraft  in  ground  effect 
using  the  methodology  described  in  Column  4.  A numerical  solution  procedure 
applied  to  considerably  simplified  forms  of  the  viscous  flow  Navier-Stokes 
equations,  has  also  been  initiated  and  belongs  to  Column  5.  Specifically, 

Bower  and  Kotansky  (ref.  17)  report  results  for  an  incompressible,  two- 
dimensional  planar  analysis  for  a rectangular  jet  issued  from  beneath  a 
plane  surface  into  a cross-flow.  Extension  to  a two-equation  turbulence 
model  and  compressible  flow  was  reported  by  Bower  (ref.  18)  and  the  discussion 
includes  the  theoretical  extension  to  three-dimensions  using  the  streamfunction 
vorticity  equation  system. 

An  alternative  approach,  also  belonging  to  Column  5 of  Table  1,  is 
developed  herein.  The  analysis  of  the  experimental  data  clearly  indicates 
that  three-dimensionality  is  absolutely  intrinsic  to  the  flowfield  develop- 
ment; hence,  there  is  no  lower  dimensional  description  that  can  yield 
meaningful  results.  The  jet  and  near-field  flow  is  rotational,  viscous  and 
undoubtedly  highly  turbulent.  Therefore,  for  a prediction  methodology  to  be 
adequate,  it  must  be  obtained  using  some  form  of  the  three-dimensional  time- 
averaged  Navier-Stokes  equations.  Furthermore,  since  the  three-dimensional 
analysis  is  required,  the  complete  Reynolds  stress  tensor  must  be  considered 
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Prediction  of  Lift-jet/ai rframe/ground  Interactions  (Ref.  15) 
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as  applicable.  Remote  from  the  plate  there  is  indication  of  an  inviscid 
flow  streamline  connecting  upstream  and  downstream,  transverse  plane  velocity 
stagnation  points.  Both  the  initial  jet  and  onset  cross-flow  velocities  are 
steady  and  the  dominant  scalar  component  of  the  jet  velocity  vector  is  uni- 
directional, parallel  to  the  curvilinear  jet  axis. 

These  observations  guided  the  analysis  conducted  under  this  contract 
which  has  yielded  derivation  of  a three-dimensional,  viscous-inviscid  inter- 
action algorithm  for  detailed  prediction  of  the  near-field  flow  associated 
with  a V/STOL  jet  injected  transverse  to  a subsonic  crossflow.  Following 
establishment  of  the  differential  equation  systems  required  solved  (in  the 
next  section)  a complete  description  of  the  sequential  steps  of  the  iterative 
interaction  algorithm  is  presented.  Table  2 provides  a brief  summary  of  the 
major  steps.  Following  Input,  a sequence  of  inexpensive  two-dimensional 
potential  flow  solutions  are  conducted  at  the  injector  surface,  which  yield 
determination  of  the  interface  contour  separating  the  potential  and  viscous- 
turbulent  solution  domains,  and  the  initial  distribution  of  entrainment 
velocity  into  the  jet  region.  An  interpolation,  guided  by  any  available  data, 
yields  the  required  three-dimensional  potential  flow  solution  domain  shape. 

Step  C.  A relatively  inexpensive  three-dimensional  potential  solution 
determines  the  first  estimate  of  the  inviscid  pressure  distribution  on  the 
plate  and  on  the  viscous  jet  region  domain  boundaries.  A modestly  inexpensive 
computational  sequence  then  yields  initialization  of  all  dependent  variables 
of  th;  three-dimensional  turbulent  flow  solution,  plus  distribution  of  the 
pressure  fields  which  inter-connect  the  two  solution  domains.  Step  E is  a 
relatively  more  expensive  solution  of  the  complete  three-dimensional  viscous 
and  turbulent  jet  flowfield  equations.  Convergence  of  the  sequence  occurs 
when  the  potential  flow  entrainment  (which  is  negative)  at  the  upper  extremity 
of  the  Step  C solution  domain  and  parallel  to  the  jet  direction,  is  in  agree- 
ment with  that  predicted  by  the  viscous  jet  region  solution,  as  affected  by 
entrainment  everywhere  perpendicular  to  the  jet  axis. 

Hence,  the  developed  interaction  algorithm  is  based  upon  existence  of 
entrainment,  both  as  a primary  motive  force  as  well  as  the  check  on  solution 
compatability.  The  output  of  the  solution  algorithm.  Step  F,  includes  the 
detailed  velocity  and  turbulence  correlation  distributions  within  the  jet 
region  and  the  induced  pressure  distribution  on  the  injector  plate.  The  jet 
region  solution  will  also  yield  a computational  estimate  of  the  entrained  mass 
flux.  The  following  section  develops  the  interaction  algorithm,  which  requires 
considerable  exposition  of  differential  equation  systems.  It  is  hoped  the 
development  is  adequate  without  being  burdensome.  The  geometric  flexibility 
of  finite  element  numerical  solution  methodology  is  employed  to  establish  the 
numerics,  as  well  as  to  generate  the  computational  results  discussed  in  the 
subsequent  section.  An  Appendix  is  included  to  record  a more  complete 
development  of  the  viscous  and  turbulent  differential  equation  systems  required 
for  the  jet  region  analysis. 
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Table  2 


SUMMARY  OF  THREE-DIMENSIONAL  V/STOL  JET  INTERACTION  ALGORITHM 


PHASE 

OPERATION 

DESCRIPTION 

A. 

Input 

Problem  specified  by: 

U*,  = crossflow  velocity 

U.,d.,r.  = jet  velocity,  diameter,  axis  distribution 

J J J 

B. 

2D  Potential 

Flow 

Computati on 

1.  Find  Cp  for  blockage  contours  on  injector  plate 

2.  Extract  entrainment  velocities  at  plate. 

C. 

3D  Potential 

Flow 

Computation 

1.  Interpolate  interaction  surface 

2.  Distribute  entrainment 

3.  Compute  inviscid  flow  plate  Cp 

D. 

2D  Viscous 

and  Turbulent 

Jet  Region 

Solutions 

1.  Initialize  jet  region  axial  velocity 

2.  Compute  transverse  plane  velocity  distributions 

3.  Compute  turbulence  correlation  distributions 

4.  Compute  interior  pressure  distribution 

E. 

3D  Turbulent 

Jet  Region 
Computation 

1.  Integrate  complete  3D  jet  region  equations  using 
Step  D.  results  for  start-up.  Yields  complete  3D 
turbulent  jet  velocity  distribution 

2.  Stop  at  top  of  3D  domain  of  Step  C. 

3.  Compare  turbulent  jet  axial  velocity  to  (negative) 
entrainment  estimated  at  Step  C.2. 

4.  Repeat  analysis  starting  at  Step  C.l.  if  required. 

F. 

Output 

1.  Complete  3D  turbulent  jet  region  flowfield. 

2.  Cp  distribution  on  injector  plate. 

3.  Prediction  of  entrainment  volume  into  jet. 
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THREE-DIMENSIONAL  INTERACTION  ALGORITHM 


Overview 

The  requirement  is  to  establish  a three-dimensional  model  of  the  V/STOL 
jet  problem,  and  the  identification  of  methodology  to  yield  a numerical 
solution.  The  only  three-dimensional  flowfield  equation  system  solvable  on 
present  digital  computer  hardware  systems  is  potential  flow.  Yet  the  potential 
theory-based  algorithms,  as  briefly  discussed  in  the  Introduction,  have  fallen 
far  short  of  adequately  describing  the  near-jet  flowTield,  primarily  due  to 
lack  of  an  adequate  jet-wake  model  (c.f.,  ref  19).  An  alternative  model  is 
developed  herein  that  relies  on  three-dimensional  potential  theory  only  to 
provide  pressure  boundary  conditions  for  a completely  viscous,  rotational  and 
turbulent  description  of  the  V/STOL  jet  and  the  near-reqion  flow.  A para- 
bolized form  of  the  steady  three-dimensional  time-averaoed  Navier-Stokes 
equations  is  derived  to  describe  this  jet  region  flowfield.  An  algorithm  is 
also  developed  to  couple  the  two  separate  differential  equation  descriptions 
within  an  interaction  formalism. 


Interaction  Surface  Definition 

The  viability  of  an  interaction  algorithm,  for  an  aerodynamics  description, 
is  based  solely  upon  identification  of  a surface  within  the  flowfield  upon 
which  inviscid  potential  flow  boundary  conditions  can  be  specified.  The 
classical  approach  for  adding  boundary  layer  effects  to  an  airfoil  analysis, 
for  example,  has  been  to  assume  the  inviscid  contour  corresponds  to  the 
physical  thickness  everywhere  augmented  by  the  displacement  surface  5* 
computed  by  the  boundary  layer  solution.  In  this  case,  6*  is  assumed  an 
inviscid  streamline,  hence  impervious  to  flow. 

This  appears  a valid  concept  for  the  V/STOL  jet  problem  only  at  sufficient 
distance  from  the  surface  of  injection.  The  available  data  appear  to  confirm 
that,  for  a circular  jet  at  A = 8 for  example,  that  this  inviscidly-impervious 
surface  develops  by  about  five  diameters  downstream  of  injection.  Figures  7 
and  8 illustrated  this  observation,  the  crucial  test  being  existence  of  an 
apparent  streamline  connecting  two  stagnation  points.  From  Figure  8,  the 
stagnation  points  exist  on  the  symmetry  plane  at  about  +4  jet  diameters  from 
the  defined  center  of  the  jet  flow  region,  and  the  region  extends  laterally 
about  the  same  distance.  Hence,  for  this  case  at  least,  the  region  of  viscous 
dominated  rotational  flow  appears  bounded  by  a nominal  circle  of  radius  equal 
to  approximately  4 jet  diameters.  Experiment  (c.f.,  ref.  11 ) indicates  that 
the  size  of  this  region  increases  as  the  jet  travels  downstream,  and  that  it 
tends  to  become  nominally  ellipical  in  cross-sectional  shape  with  major  axis 
perpendicular  to  both  the  cross-flow  velocity  vector  and  the  normal  to  the 
plate. 


Adequate  vortex  models  exist  to  describe  this  far-field  evolution,  roll-up, 
and  turning  parallel  to  the  cross-flow,  see  for  example  reference  12.  The 
present  requirement  is  to  describe  the  flow  in  the  first  (say)  4-5  jet  diameters 
downstream  of  injection.  The  sole  alternative  guidance  available  for  definition 
of  an  interaction  surface  is  the  measured  surface  Cp  and  oil  streakline  photo- 
graphs. Recalling  Figures  3-4  and  the  pertinent  discussion,  no  inviscid  flow 
streamline  is  indicated.  In  particular,  mass  is  specifically  injected  into  the 
jet  wake  from  the  cross-flow.  Fortunately,  and  of  specific  usefulness  here, 
the  interaction  surface  separating  the  potential  and  viscous  flow  domains  need 
not  be  an  inviscid  streamline. 

The  differential  equation  governing  potential  flow  is  simply  the  linear 
Laplacian 


L(<j>)  i v2<f>  = 0 


where  the  dependent  variable  is  the  (perturbation)  potential  function  <f>, 
defined  in  term^  of  the  local  velocity  vector  u,  and  the  reference  (cross- 
flow)  velocity  as 


u = U„  - v*  (2) 

Since  equation  (1)  is  an  illiptic  boundary  value  problem,  boundary  conditions 
on  <t>  are  everywhere  required  specified.  If  the  boundary  corresponds  to  an 
inviscid  streamline,  than  u must  there  lie  tangential.  The  equivalent 
expression  is  that  the  dot  product  of  u with  the  normal  vector  n to  the  surface 
vanishes,  i.e.  u.n  = 0.  From  equation  (2),  the  computational  equivalent  is 

v<t>.n  = iL.n  (3) 


which  states  teh  normal  projection  of  equals  the  onset  component  of  U^,. 

Should  the  boundary  surface  not  be  an  inviscid  streamline,  then  ij.h  does 
not  equal  zero  but  instead  equals  the  inviscid  "onset"  velocity  component  UL  . 
This  is  the  component  of  Tj,  equation  (2),  that  passes  through  the  surface. 
Therefore,  the  boundary  condition  for  solution  of  the  potential  distribution 
<p  becomes 


v<j> . n = Uco-h  - ul 


(4) 


Use  of  equation  (4),  which  is  simply  a modest  generalization  of  equation 
(3),  as  the  boundary  specification  for  solution  of  equation  (1)  for  potential 
flow,  yields  the  required  flexability  for  the  V/STOL  interaction  algorithm. 

For  the  streaklines  shown  in  Figure  3a,  corresponding  to  solid  blockage,  the 
visible  inviscid  streamline  not  surprisingly  appears  in  the  nominal  shape  of 
an  airfoil.  Therefore,  assume  the  blockage  effect  for  the  V/STOL  jet, 

Figure  3b,  is  also  an  airfoil  contour,  the  thickness  ratio  and  chord  of  which 
is  to  be  determined.  It  probably  is  not  totally  different  from  that  of  the 
solid  blockage  contour,  but  it  must  be  determined.  The  required  data  for 
this  is  specifically  available  from  the  experimental  measurements  of  surface 
Cp  distributions  for  many  configurations. 

A planar  airfoil  contour  Cp  is  readily  computed  for  no  "onset"  flow, 
i.e.  no  entrainment  using  equations  (l)-(3).  Since  the  real  flow  with 
entrainment  must  uniformally  involve  equal  or  larger  velocities,  for  any 
selected  airfoil  shape,  the  experimental  Cp  must  be  either  equal  to  or 
more  negative  than  the  zero  entrainment  case.  The  minintm  tMckness  (ratio) 
for  the  no  entrainment  case  is  the  jet  diameter,  but  no  upper-bound  exists 
on  the  thickness.  Therefore,  it  is  reasonable  to  expect  that  some  blockage 
ccntour(s)  can  be  found  nreeting  these  requirements.  Assuming  this  valid 
(as  is  indeed  verified  in  the  next  section  for  the  round  jet  case  with  x = 8), 
'■.he  definition  of  pressure  coefficient  Cp  yields  the  local  distribution  of 
entrainment  velocity  u1-  as 


C 


P 


(o1-)2 


(5) 


where  C^  is  the  computed  blockage-only  distribution.  The  direction  of  u1- 

is  of  course  inwards;  the  availability  of  n allows  resolution  of  u1-  into 
appropriate  scalar  components  of  the  local  velocity  vector  u. 


This  procedure  yields  definition  of  a (range  of)  inviscid-viscous  flow 
interaction  contour(s)  on  the  injection  plate  surface.  The  shape  of  the 
blockage  surface  is  also  known  (for  the  round  jet  with  x = 8,  at  least)  at 
some  distance  from  the  plate.  Additional  experimental  data  would  allow 
definition  of  the  interaction  surface  between  these  extrema.  In  its  absence, 
it  is  fair  to  assume  that  the  surface  would  smoothly  cornect  the  two  curves 
and  would  have  to  be  determined  by  the  analysis.  This  requires  that  an 
iterative  procedure  be  developed  between  the  potential  and  viscous  flow 
solutions  with  interface  shape  and  entrainment  the  iteration  parameters. 

This  leads  to  derivation  of  the  key  feature  of  pressure  coupling  within  the 
interaction  algorithm. 
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Pressure  Coupling  Algorithm 


A parabolized  form  of  the  steady,  three-dimensional  time-averaged  Navier- 
Stokes  equations  will  be  employed  for  the  jet-wake  region  flowfield  as  discussed 
in  the  following  section.  Bearing  this  in  mind,  a key  element  in  the  paraboli- 
zation procedure  is  that  the  transverse  plane  pressure  distribution  and  axial 
pressure  gradient  be  locally  solvable  as  a function  for  distance  along  the  jet 
trajectory.  This  is  accomplished  by  derivation  of  a differential  equation 
obtained  from  the  divergence  of  the  transverse  plane  mean  flow  momentum 
equations.  The  detailed  development  is  presented  in  the  Appendix;  it  is 
sufficient  at  this  juncture  to  simply  state  the  derived  differential  equation. 
Defining  the  limited  tensor  index  summation  convention,  1 < j <_  3 and 
2 < k,  i <_  3,  where  x£  is  the  coordinate  system  spanning  the  plane  transverse 
to  the  local  jet  direction,  with  unit  normal  parallel  to  the  xx  coordinate, 
the  differential  equation  governing  the  jet-region  pressure  distribution  is 


Lfp>  ■ w 


S£_ 

3x„ 


nr 

J 


3^  (5£k  ' PU£Uk} 


= 0 


(6) 


The  overbar  and  over-tilde  in  equation_(6)  signify  time-averaged  and  mass- 
weight  time-average  respectively,  and  p is  pressure,  p is  density,  u^  is 
velocity  vector,  is  the  Stokes  viscous  stress  tensor,  and  pu£u' 

is  the  Reynolds  stress  tensor. 

The  solution  algorithm  for  equation  (6)  is  based  upon  the  observation  that 
the  solution  of  a linear,  elliptic  Poisson  equation  consists  of  a complementary 
and  a particular  contribution  as 


P(xi)  = pc(xi)  + Pp(x.)  (7) 

By  definition,  the  complementary  solution  satisfies  the  homogeneous  form  of 
equation  (7) , i .e. , 


l(Pc)  = 


3X2 


= 0 


(8) 


subject  to  the  known  boundary  conditions  for  pressure  p(x-).  Since  the  inter 
face  invisc id  flow  pressure  distribution  is  everywhere  kn&wn. 


In  equation  (9),  x^  indicates  constrained  to  the  interfacial  boundary  of 
the  viscous  solution  domain,  and  pi  is  the  (inviscid  flow)  pressure  level 
distribution.  Elsewhere,  on  the  jet  region  central  symmetry  plane  specifically, 
the  appropriate  boundary  condition  is  vanishing  normal  gradient. 


*(PC)  = 


(10) 


Therefore,  the  solution  to  equations  (8)-(10)  provides  determination  of  the 
pressure  distribution  interior  to  the  viscous  jet  region  induced  by  the 
impressed  cross-flow  velocity  (distribution).  Since,  as  presented  in  the 
previous  section,  the  inviscid  flow  pressure  distribution  is  a priori  known, 
the  solution  to  equation  (8)  provides  the  required  pressure  coupling  of  the 
interaction  algorithm. 

The  particular  pressure  distribution  pp  modifies  the  computed  complemen- 
tary pressure  distribution  to  account  for  viscous  modifications  within  the  jet 
region.  The  particular  solution  pp  is  any  function  satisfying  equation  (6) 
with  homogeneous  boundary  conditions  on  interface  segments  coincident  with 
pi  known.  Since  the  convection  term  in  equation  (6)  vanishes  identically  on 
a symmetry  plane,  the  appropriate  form  of  the  boundary  condition  is 


*(Pp)  = 


ap  - 
3xjt 


puku; 
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Since  the  second  term  in  equation  (11)  obliterates  the  boundary  condition 
stress  term  in  equation  (6),  the  appropriate  Poisson  equation  for  pp  must  be 
the  inviscid  form. 
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puju* 


pVk 


= 0 


(12) 


Equations  (6 ) - ( 12 ) define  the  algorithm  for  determination  of  pressure 
distributions  in  the  three-dimensional  parabolic  viscous  flow  solution.  The 
pressure  solutions  are  summed,  as  in  the  defining  equation  (7)  and  differen- 
tiated to  form  the  transverse  pressure  gradient  field.  This  completes  defini- 
tion of  the  pressure  coupling  algorithm  for  the  interaction  algorithm,  and 
attention  is  now  focused  on  establishment  of  the  viscous  jet  region  differential 
equation  description. 
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Jet  Region  Governing  Differential  Equation  System 

The  singular  requirement  for  the  proposed  viscous-invisid  interaction 
algorithm  and  a significant  contribution  of  this  study  is  derivation  and 
development  of  a numerically  tractible  differential  equation  system  describing 
the  three-dimensional  viscous  and  turbulent  flowfield  within  the  jet  region 
interior  to  the  interaction  interface.  The  developed  differential  equation 
system  is  named  the  three-dimensional  parabolic  Navier-Stokes  (3DPNS) 
equations,  and  has  emerged  after  several  years  development  in  a consequentially 
complete  form.  The  3DPNS  system  is  obtained  from  the  complete,  non-linear 
time-averaged  Navier-Stokes  partial  differential  equation  set,  following 
application  of  simplifying  assumptions  which  are: 

(1)  The  flowfield  is  steady 

(2)  A dominant  direction  for  the  mean  flow  is  uniformly  discernible 

(3)  In  this  direction  (only)  diffusion  is  negligible  compared  to  all 
other  processes  including  convection,  and 

(4)  The  overall  three-dimensional  elliptic  boundary  value  character 

is  enforced  by  interaction  with  a concomitantly  computed  potential 
flow 

The  examiniation  of  the  character  of  the  V/STOL  jet-region  flow,  as 
discussed  to  the  present,  indicates  that  each  of  these  assumptions  is  satisfied. 
The  flowfield  can  be  assumed  steady,  and  because  of  existence  of  the  plate  of 
injection,  the  jet  and  jet-wake  flows  are  uniformly  directed  therefrom  along 
the  jet  axis.  (The  sole  exception  appears  to  occur  at  the  rear  stagnation 
region,  as  evidenced  in  the  oil  streakline  photos,  see  Figure  3b  and  4.  The 
jet-region  solution  domain  therefore  can  at  most  intersect  this  location,  but 
cannot  include  it.  This  places  a limitation  on  the  analysis  that  can  be 
removed  only  by  going  to  a complete  three-dimensional  Navier-Stokes  descrip- 
tion which  is  computationally  intractible  at  present.)  For  sufficiently  large 
jet  velocity  ratio  the  third  assumption  appears  valid,  and  the  entire 
development  is  derived  within  an  interaction  algorithm  with  inviscid  exterior 
three-dimensional  flowfield. 

Establ ishement  of  the  3DPNS  differential  equation  system  is  a fairly 
consuming  process  and  the  derivation  is  presented  in  completeness  by  Baker 
(ref.  20  ).  The  Appendix  of  this  report  presents  a summary  of  the  developmental 
procedures.  For  the  purpose  of  completing  the  interaction  solution  algorithm 
definition,  the  3DPNS  differential  equation  system  that  governs  convection, 
diffusion,  generation  and  annihilation  of  both  mean  flow  and  turbulence 
correlation  phenomena  within  the  jet-wake  solution  domain  lying  interior  to  the 
interaction  surface,  is 
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Recall  the  index  limited  summation  convention,  1 < i <3  and  2 s k,  1 s 3, 
and  that  the  variables  appearing  in  equations  ( 1 3 ) - ( 1 8 ) are  the  mass-weighted, 
time-averaged  density,  pressure  and  velocity  and  dynamic  Reynolds  stress  tensor. 
Re  is  the  reference  Reynolds  number.  Equations  (17)- (18)  are  the  three- 
dimensional  parabolic  form  of  the  transport  equations  for  turbulence  kinetic 
energy  k and  isotropic  dissipation  function  e,  defined  as 


k 


(19) 


2 . _ 3u[ 

3 6UG  = 


(20) 


These  equations  are  required  solved,  in  concert  with  equations  (13)-(16),  to 
complete  closure  of  the  3DPNS  system  for  the  Reynolds  stress  distribution. 

The  dynamic  Reynolds  stress  is  assumed  the  product  of  the  time-averaged  density 
and  the  kinematic  Reynolds  stress  tensor,  the  parabolic  components  of  which  in 
principal  coordinates  are  (ref.  21) 
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Finally,  the  various  C-j  and  C“  are  correlation  constants  of  the  turbulent 
closure  model,  see  the  Appendix. 


Equations  ( 14)- ( 18)  are  five  initial-value,  two-dimensional  elliptic 
boundary  value,  non-linear  partial  differential  equations  written  on  the 
dependent  variables  of  jet  region  velocity  u-j  and  turbulence  correlations  k and 
e.  Equation  (13)  is  a constraint  equation  enforcing  conservation  of  mass, 
which  is  inserted  directly  into  the  right-hand  non-homogeniety  of  the  particu- 
lar pressure  Poisson  equation  (12).  The  axial  pressure  gradient  distribution 
3p/3xj , used  to  drive  the  dominant  jet-flow  momentum  equation  (14),  is  estab- 
lished solely  from  the  complementary  pressure  solution.  The  total  computational 
pressure  field  is  employed  to  drive  the  transverse  plane  velocity  components, 
equations  ( 15) - ( 16 ) . 

The  solution  domain  for  equations  ( 14)- ( 18)  is  the  three-dimensional 
region  interior  to  the  interface  surface  and  closed  from  below  by  the  injection 
plate,  see  Figure  10.  The  domain  is  open-ended  in  the  downstream  jet  direction, 
the  result  of  the  parabolizing  simplification.  Hence,  the  3DPNS  system  is 
initial-value  on  the  xx  curvilinear  coordinate,  Figure  10,  and  the  solution  is 
started  at  or  near  the  injection  plate  surface.  As  a consequence,  initial- 
condition  distributions  are  required  specified  for  all  dependent  variables 
(q > = {p,  Dj,  k,  e}  at  this  location.  This  in  itself  is  a consequential  compu- 
tational problem  since  in  reality  all  that  is  known  is  the  (nominal)  jet  velocity 
Uj,  on  a circle  of  diameter  d,  and  the  distribution  of  the  tangential  and 
entrainment  velocity  components  of  u A on  the  interface  surface.  Procedures  have 
been  developed  to  meet  the  differential  equation  requirements,  however,  and 
are  discussed  with  the  numerical  results.  Boundary  conditions  are  required 
specified  everywhere  on  the  (two-dimensional)interface  surface  to  render 
equations  ( 14>(  18)  a well-posed  problem  description. 
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At  the  solution  initiation  plane,  and  as  a consequence  of  the  iterative  inter- 
action algorithm,  the  distribution  of  the  tangential  and  entrainment 
components  of  the  transverse  plane  velocity  field  u^  are  everywhere  known  from 
the  exterior  potential  flow  Cp  distribution.  The  boundary  conditions  for  ui 
can  be  similarly  extracted  from  the  exterior  potential  solution  since  a non- 
zero axial  (Xl)  pressure  gradient  induces  a corresponding  parallel  velocity 
component.  The  boundary  conditions  for  k and  e are  specified  either  to  agree 
with  the  freestream  turbulence  level,  if  known,  or  as  a vanishing  normal  deriv- 
ative. Although  the  theory  is  conceptially  extendible  to  compressible  and/or 
heated,  dissimilar  jets,  the  current  analysis  is  restricted  to  constant  uniform 
density  flows. 


Coordinate  System  Descriptions 

Experimental  evidence  confirms  that  the  cross-sectional  area  and  shape  of 
the  V/STOL  jet  interaction  surface  is  strongly  a function  of  distance  along 
the  jet  path.  For  the  circular  jet  at  X = 8,  for  example,  the  cross-section 
changes  from  an  approximate  35%  thick  airfoil  at  the  injection  surface  to 
almost  circular  at  xj/d  = 5.3.  A requirement  of  the  3DPNS  analysis  of  the 
jet-region  flow  evolution,  therefore,  is  to  match  a generally  non-regular 
shaped  solution  domain  closure  that  changes  with  solution  progression  down- 
stream. Furthermore,  since  xx  is  a curve,  the  entire  problem  must  be  cast  in 
an  appropriate  curvilinear  coordinate  system. 

The  overall  requirements  are  met  by  employing  a finite  element-based 
numerical  solution  algorithm  for  solving  the  3DPNS  equations.  Even  with  this 
inherent  geometric  flexibility,  a grid  stretching  coordinate  transformation  is 
required  employed  such  that  the  transverse  plane  discretization  continuously 
spans  the  entire  interior  bounded  by  the  interaction  surface.  Figure  10 
illustrates  the  global  x,  y,  z coordinate  system,  with  origin  at  the  initial 
jet  center,  and  the  curvilinear  x^  system.  The  xj  coordinate  is  defined  as 
everywhere  parallel  to  the  jet  centerline,  defined  as  the  locus  of  extremum 
Ui  lying  on  the  symmetry  centerplane  (c.f.,  ref.  12).  For  simplicity  of 
illustration,  assume  that  the  jet  region  interface  surface  is  required  to 
stretch  only  in  the  lateral  direction,  i.e.  parallel  to  x2.  A coordinate 
stretching  transformation  that  will  continuously  place  the  lateral  span  between 
the  jet  symmetry  plane  and  the  interface  surface  onto  a constant  length  is 

x2  - f i ( x i ) 

n = 

[f2(xi)  - fiUojf"1  (22) 

In  equation  (22),  fi(xi)  is  the  symmetry  plane,  f2(xi)  is  the  maximum  lateral 
span  of  the  interface  surface  as  a function  of  location  Xl,  and  f is  a nor- 
malizing factor.  The  chain  rule  for  differentiation  yields 
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where  superscript  prime  denotes  the  ordinary  derivative  with  respect  to  xl 
Hence,  in  terms  of  coefficients  h-j . 1 <_  i £ 3, 


fif2  - f,)-1 

= - f1h1 

f-Hf'z  - f>i 


related  to  the  metric  of  equation  (22),  the  spatial  derivatives  in  3DPNS  are 
reexpressed  as 


3Xj  35 


[h2  + nh3j: 


3 . h 3_ 

3x2  ~ nj3n 


Then  computations  would  be  performed  in  the  fixed  5,  n,  x3  coordinate 
system. 

In  general,  the  interface  surface  will  expand  in  both  the  x2  and  x3 
coordinate  directions  as  the  jet  proceeds  downstream.  The  required  generali- 
zation of  equation  (22)  is 


*2  “ ^21 ^ X1 ^ 
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x3  " f31*xl* 


In  equation  (26),  f2i  are  the  counterparts  of  f-j  in  the  one-dimensional  case 

while  f3-j  describe  the  extremum  span  of  the  interface  surface  parallel  to 

the  x3  coordinate  direction.  Again  employing  the  chain  rule  for  differentia- 
tion on  x-j  yields 
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The  functions  h^,  1 £ i <_  3,  are  again  related  to  the  metric  of  the  coordinate 
transformation,  equation  (26) , and  are  defined  as 
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The  superscript  prime  denotes  the  ordinary  derivative  with  respect  to  Xj. 

The  nj  coordinate  system  is  fixed  in  the  transform  space  and  hence  insensitive 

tO  X]  . 


One  additional  transformation  is  required  to  account  for  Xj  being  curvi- 
linear. Since  in  the  computational  solution,  integration  of  the  3DPNS  equations 
occurs  step-by-step,  displacements  parallel  to  the  xt  coordinate  are  re-expressed 
in  terms  of  an  angular  displacement  e as 


dxj  = r^Jde 


(29) 


In  equation  (29),  r is  the  local  radius  to  the  jet  centerline,  and  0 is  measured 
from  the  negative  x axis  of  the  global  coordinate  system,  see  Figure  10.  Ample 
experimental  data  are  available  to  evaluate  r(xj)  for  circular,  square  and 
rectangular  V/STOL  jets  over  a range  of  velocity  ratios.  For  example,  Figure  11a 
illustrates  the  curves  determined  by  McMahon  and  Mosher  (ref.  6)  for  circular  and 
rounded-corner  rectangular  jets.  Figure  lib,  from  reference  12,  illustrates  the 
collapse  of  circular  jet  data  onto  a universal  curve  obtained  using  the  jet 
velocity  ratio  as  a similarity  parameter.  Kamotani  and  Greber  (ref.  11)  deter- 
mined that,  using  their  data  for  the  circular  jet  for  4 £ x £ 10,  the  equation  of 
the  jet  centerline  for  all  cases  was  expressible  as 
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a)  Various  Jet  Shapes  (Ref.  6) 


b)  Circular  Jet  at  Various  A (Ref.  12) 


Figure  11.  Measured  Jet  Centerline  Trajectories 
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The  x and  z are  global  coordinates  and  d is  the  diameter  of  the  jet  orifice 
in  equation  (30) 

The  3DPNS  equation  system  ( 13)- ( 18)  must  bo  expanded  into  the 
x-j  e { x i , x2,  x3}  e (e,  z,  r}  coordinate  system,  prior  to  employing  the  grid 
stretching  transformations  described.  The  details  are  messy,  but  the  end 
point  modifications  are  modest.  For  example,  the  xj  (dominant  flow)  momentum 
equation  (14)  becomes 
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There  is  an  added  convection  term  and  the  total  stress  tensor  t.  has  the 
scalar  components  in  {e,  z,  r}. 
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The  grid  stretching  transformation  parallel  to  the  x2  coordinate  is  unaffected. 
Viewing  equation  (22)  for  simplicity,  however,  and  noting  that  x3  = r in  the 
present  case,  the  transformation  is 

r - f ! (e ) 

’ LMeKjMf1  (33) 


The  f-j(o)  describe  the  transverse  extent  of  the  interface  surface  parallel  to 
the  r axis  and  are  assumed  expressed  as 
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discussed,  and  noting  that  the  renuired  derivative  always 
’3/39,  the  polar  coordinate  system  eouivalent  of  eauations  (25) 
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The  hj  defined  in  equation  (35)  are  identical  to  those  in  equation  (24)  but 
divided  by  r.  In  actual  practice  the  inverse  r modifying  3/36'. can  be 
replaced  by  the  constant  r for  sufficiently  large  r,  as  occurs  in  the  near- 
field jet  region,  for  example. 


Interaction  Solution  Algorithm 

The  exterior  potential  flow  and  interior  3DPNS  differential  equation 
systems,  required  solved  for  the  V/STOL  jet  interaction  algorithm,  have  been 
identified.  The  boundary  separating  these  two  solution  domains  is  the  inter- 
action surface  upon  which  computed  flowfield  variables  must  come  into  agree- 
ment. The  geometry  of  the  intersection  of  the  interaction  surface  with  the 
injection  plate  is  assumed  known,  as  determined  by  comparison  of  blockage  plus 
entrainment  pressure  predictions  with  experimental  data.  The  radii  of  the 
family  of  circles  locally  tangent  to  the  jet  trajectory  is  also  known  from 
experimental  data.  In  many  instances,  for  example,  the  circular  jet  at  X = 8, 
the  geometry  of  the  intersection  of  the  interaction  surface  with  the  trans- 
verse plane,  perpendicular  to  the  jet  axis,  is  also  known  at  one  or  more  down- 
stream stations.  The  interaction  surface  geometry  between  the  injection  plate 
and  downstream  intersections  is  typically  unknown. 

Even  if  the  interaction  surface  description  were  everywhere  known,  an 
iteration  procedure  is  required  developed  between  the  3DPNS  and  3D<f>  solutions. 
The  basic  iteration  parameter  for  assessment  of  convergence  is  the  accelera- 
tion (parallel  to  the  jet  axis)  of  the  fluid  contained  within  and  entrained 
into  the  3DPNS  solution  domain.  Figure  12  illustrates  the  symmetric  half- 
domain employed  for  the  3D<j>  solutions  discussed  in  the  next  section.  The  3D<t> 
solution  requires  the  onset  velocity  be  everywhere  specified  on  the  domain 
boundary,  recall  equations  (l)-(4).  This  includes  the  top  or  "lid"  of  the 
domain,  shown  as  a plane  parallel  to  the  global  x,  y plane.  The  loss  of  mass 
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to  the  3DPNS  solution  domain  is  termed  the  entrainment  u-1-  . As  a consequence, 
an  acceleration  of  Uj  is  distributed  throughtout  the  3DPNS  interior  by  the 
action  of  continuity,  equation  (13).  Hence,  the  3DPNS  solution  will  yield 
u j >0  selectively  at  the  3D*  "lid"  intersection  with  the  3DPNS  domain,  see 
Figure  12.  This  results  in  a non-zero  "onset"  velocity  boundary  condition 
for  3D*,  which  can  be  distributed  across  the  lid  by  interpolation  to  zero  at 
the  lateral  extremity. 

Hence,  the  action  of  entrainment  from  the  3D*  solution  into  the  3DPNS 
solution  becomes  reflected  c's  an  onset  velocity  boundary  condition  modification 
for  3D*.  The  two  solutions  are  accepted  when  this  difference  in  onset  velocity 
components  vanishes  to  within  a specified  tolerance.  Since  the  elevation  of 
the  3D*  lid  above  the  plate  is  not  required  to  be  farfield,  the  interaction 
solution  iteration  is  presumably  operable  on  progressively  larger  domains  as 
generated  by  addition  of  layers  to  the  3D*  domain.  The  iteration  sequence  is 
presented  in  complete  detail  in  Table  3.  The  previous  discussion  on  the  con- 
vergence test  occurs  at  step  E.2.  Note  that  an  additional  iterative  decision 
is  suggested  at  step  C.3  prior  to  actual  generation  of  a 3DPNS  solution.  The 
sole  input  parameters  are  freestream  and  jet  initial  velocity,  jet  diameter, 
and  radius  distribution  of  the  jet  centerline.  All  remaining  steps  require 
completion  of  various  computational  sequences.  Therefore,  appropriate  numeri- 
cal solution  algorithms  are  required  identified. 


Finite  Element  Solution  Algorithm 


The  3DPNS  partial  differential  equation  system,  oroposed  to  describe  the 
viscous  and  turbulent  flowfield  in  the  immediate  vicinity  of  a V/STOL  jet  in 
crossflow,  has  been  presented.  Each  member  of  the  equation  set  (13)-(18) 
plus  (8)  and  (12  ) for  the  complementary  and  particular  pressure  fields, 
belongs  to  the  general  class  of  second-order,  non-linear  elliptic  oartial 
differential  equations.  Hence,  define  the  generalized  dependent  variable 
descriptor  {q > = , k,  e,  pc,  pp,  *},  where  * is  the  potential  function 

describing  the  exterior  three-dimensional  potential  flow  which  comnletes  the 
interaction  algorithm.  Each  governing  partial  differential  equation  for 
(q } belongs  to  the  form 
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“l*  3*7 

(36) 


In  equation  (36),  the  tensor  indices  range  2 <_  z <_  3,  and  1 _<  i.  _<  3,  K is  the 
diffusion  coefficient,  f!  is  a function  of  its  argument  that  specifically 
includes  three-dimensional  convection,  p is  any  solution  parameter  including 
another  dependent  variable,  and  f2  is  the  initial-value  operator  if  present. 
The  three-dimensional  solution  of  each  3DPNS  member  of  {q}  is  required 
established  on  the  bounded  open  domain  q,  formed  by  the  interior  to  the  inter- 
action surface  truncated  by  the  injector  plate  and  open  at  the  downstream  end. 


Table  3 


ITERATION  SEQUENCE  FOR  V/STOL  INTERACTION  ALGORITHM 


A.  SOLUTION  PARAMETERS 


Input: 

Uco 

= freestream  velocity 

Uj 

= initial  jet  velocity 

d 

= initial  jet  diameter 

r(x{, 

) = distribution  of  jet  axis  radius 

B.  GROUND  PLANE  CONTOUR  AND  EDGE  ENTRAINMENT 

Determine:  Si  = interaction  surface  intersection  with  plate 

u1  3 entrainment  velocity  on 

Procedure:  a)  solve  V2<f>  = 0 on  2D  blockage  "airfoil"  contours 

b)  select  Sj  with  Cp>_Cp  everywhere 
Yields : Sj  = initial  interface  contour 

u1  = initial  contour  entrainment  velocity 
Sj  = location  of  jet  inside  Sj 

. C.  3-D  POTENTIAL  FLOW 

1.  Determine:  (x x ? interaction  surface  contours  at  n stations,  starting  at 

plate  to  "sufficienctly"  downstream,  1 < i <_n . 

Procedure:  a)  using  Sj  and  any  available  data  for  Sn, 
interpolate  a smooth  family  of  contours. 

2.  Determine:  CD(S^)  = exterior  potential  flow  pressure  on  S-j(xi) 

r I 

Procedure:  a)  interpolate  u (x  ) over  S-j , 1 < i £n 

b)  apply  as  onset  velocity  boundary  conditions 

c)  estimate  ui  on  lid  of  3D  domain 

d)  solve  V2<J>  = 0 on  3D  domain 

Yields:  Cp(s-j) 

Cp  on  injector  plate 

3.  Decision:  compare  CR  on  plate  to  experimental  data (if  available) 

a)  if  agreement  good,  proceed  to  step  D 

b)  if  not,  return  to  C.l,  alter  and/or  u1 
interpolation 
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D.  3DPNS  SOLUTION  INITIATION 


1.  Determine:  u 


2.  Determine: 


3.  Determine:  p 


E.  3DPNS  JET  REGION 
1.  Determine: 


2.  Decision: 


j(xv,x")  - initial  velocity  distribution  interior  to  Sj 

Procedure:  a)  S.(x„)  = °"  si  fns1de  S> 

1 * (e > 0 elsewhere  inside  Sj 

b)  solve  V2U£ = 0 interior  to  Sj  boundary 
conditions  provided  by  resolution  of 
ub  and  u1  onto  u^x^) 

Yields:  Gi  (xji  »xi ) 

k,c  5 remaining  ini tial -value  variable  distributions  on  Sj 
Procedure:  a)  determine  VU!  • VU].  on  S^x^) 

b)  using  MLT  concepts,  determine  ve 

c)  defining  length  scale,  solve  for 
k(x£),  e(xA) 

d)  set  background  level  of  k > 0 

Yields:  k(x^,x?) 

e(xrx?) 

ve  (x^,x°) 

..(x^.x  ) = complementary  pressures  on  S. , 1 < i < n 
Procedure:  a)  solve  V2pc  = 0 on  each  S-j 

b)  store  computed  pressure  distributions 
Yields:  Mx£’xi) 

Vpc(x£,x1) 


SOLUTION 

3DPNS  solution  interior  to  S(x x ) 

Procedure : a)  integrate  3DPNS  equation  system  on  x} 

b)  all  boundary  conditions  are  known 

c)  start  pp  solution  when  derivatives  available 

d)  entrainment  distributed  by  continuity 

Yields:  ijj.u^.k.e.p  throughout  S(Xj)  envelope 

compatabi  1 ity  of  3DPNS  and  3D<t  solutions  at  Sn 
Procedure^  a)  3DPNS  yields  Uj(x^,Xj)  > 0 on  Sp 

b)  compare  this  distribution  to  potential 
lid  Uj  estimate,  C.2.c). 

Choice:  a)  if  agreement  good,  solution  is  complete 

b)  if  not,  return  to  C.2  for  next  iteration 
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The  three-dimensional  potential  flow  is  required  determined  exterior  to  the 
interaction  surface  and  above  the  injector  plate  with  farfield  closure  located 
a sufficiently  remote  distance.  Expressed  concisely,  the  3DPNS  solution  domain 
is 


S2  = R2X 


Xl) 


(37) 


while  the  3D<)>  solution  domain  is  simply  R3ex1- . 

Note  in  equations  (12)  and  (8),  for  pp  and  pG,  that  f2  vanishes  identically 
and  Xi  appears  only  as  a parameter.  Independent  of  this,  each  3DPNS  equation 
is  an  elliptic  boundary  value  description  and  boundary  conditions  are  required 
specified  on  the  interface  surface  at  each  xx  station.  As  a reminder,  the 
closure  3ft  of  the  3DPNS  solution  domain  ft  is 


3fi  = 3R2x[xq,  x)  £ y x ! (38) 

On  3n,  a boundary  condition  statement  that  includes  all  admissible  combinations 
of  function  and  normal  derivative  constraints  is 


£(q)  = a:  q + a2  K-^3-  n£  + a3  = 0 (39) 

where  the  a^Xj.Xj)  are  specified.  Except  for  q = {pc,  pp}  an  initial 
condition  is  required  at  or  near  the  injection  surface,  i.e.  on  R2U3R2  * x“. 

The  form  is 

q(x°,  x2,  x3)  e qQ(x2,  x3)  (40) 

The  finite  element  numerical  solution  algorithm  for  3DPNS  simply  assumes 
all  q and  p interpolated  on  the  discretization  of  the  solution  domain  n in 
the  form 

qe(Vx)  = {Nk(x*)}T{Q(x)}e  (4U 

where  the  subscript  e denotes  constrained  to  the  finite  element  domain 

= * xj.  Determination  of  the  unknown  expansion  coefficients  {q}  = I {Q}  , 

is  accomplished  using  the  Method  of  Weighted  Residuals  (ref.  20)  as 
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{N(x, ) } L(qe)dT  - x 


{N(x,)}«,(qp)dr 


R2 

e 


R2rt3R2 


* {0} 


(42) 


where  Se  is  the  assembly  operator.  The  assembled  finite  element  algorithm 
for  representative  partial  differential  equation  (36),  plus  boundary  condition 
statement  (39)  is 


(43) 


The  subscripts  e again  denote  functions  defined  on  the  finite  element 
discretization  ne  of  the  global  domain  n. 


The  rank  of  equation  (43)  is  identical  with  the  total  number  of  node 
points  on  R2UaR2  at  which  the  dependent  variable  q requires  solution.  For 
f 2 non-vanishing,  equation  (43)  is  a system  of  first-order  ordinary  differen- 
tial equations.  For  f?  vanishing  identically,  it  is  large  order  algebraic, 
and  the  matrix  structure  is  symmetric,  sparse  and  banded.  Hence,  standard 
finite  element  matrix  solution  methods  are  applicable  for  determination  of 
Pc,  pp,  and  <j> . 


An  implicit  numerical  integration  is  employed  to  solve  equation  (43)  for 
all  initial-value  problem  descriptions  in  3DPNS.  The  explicit  form  of 
equation  (43)  is,  in  this  instance, 


rc]efq>e  + <Me  + Me){Q}e  + {f\ 


- (0) 


(44) 


The  superscript  prime  denotes  the  ordinary  derivative  with  respect  to  the 
X!  (m ) coordinate.  The  matrix  expressions  for  the  first  three  terms, 
which  account  for  axial  convection,  and  transverse  plane  convection  and 
diffusion  respectively,  are 
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w. 

Me 

M. 


[pu,  }, 


«».>J 


(N)  f Nf  INl’dr 


{N}  {N}  dr 

3X£ 


(cve 


mm  imdT 


Dx„  3X„ 
£ £ 


(45) 


Summation  on  2 _<  i <_  3 is  implied,  and  derivatives  on  x-j  are  replaced  by 
equations  (26)  in  the  sequel.  All  terms  not  explicitly  involvinq  q,  i.e. 
f Q } , are  contained  in  {f>e  . 

A single-step,  implicit  numerical  integration  algorithm  is  employed. 


■ 

{Q}j  ♦ h|e{Q}j+1  ♦ (i  . e){Q}: 


(46) 


In  equation  (46),  j is  the  jet  axis  station  index,  h is  integration  step-size 
Axl5  and  9 equals  one-half  yields  the  trapezoidal  rule.  Following  the  usual 
matrix  manipulations,  insertion  of  equation  (44)  into  (46)  yields  a large 
order,  non-linear  algebraic  equation  system.  The  Newton  matrix  iterative 
algorithm  for  solution  of  this  system  is  employed  as 


(Q} 


M 


(47) 


The  dependent  variable  in  equation  (47)  is  the  iteration  vector,  and 


* <SQ,jt!  <«> 

where  p is  the  iteration  index,  and  i < 0 is  an  integer  that  can  retard 
evaluation  of  the  Jacobian  as  an  economy  measure.  The  right  side  of  equation 
(47)  is  the  homogeneous  form  of  equation  (46)  evaluated  with  the  pill  iterate, 
i.e. 


38 


J 


1 


(49) 


'F>S*1  • Se 


“ 

’ ' 

f 

* 

«>5.1  - {q>Ji 

* hle{ge}j+1  + (1  ' 9){ge}j) 

J 

where 


(9e>;  4 <[()],♦  [K],)  (Q)J  * (f), 


(50) 


Equations  (49)- (50)  are  defined  only  in  terms  of  inner  products  on  elements, 
with  the  assembly  operator  yielding  the  equivalent  global  expression.  The 
vanishing  of  { F } , to  within  definition  of  a computed  zero,  yields  equation 
(47)  homogeneous,  hence  convergence  of  the  iteration  for  any  evaluation  of 
the  Jacobian.  The  initial  estimate  (Q}j+1  for  any  iteration  can  be  deter- 
mined using  e = 0.  J 

By  definition,  the  Jacobian  is  the  derivative  of  equation  (49)  with 
respect  to  (Q }p.  It  is  computationally  evaluated  as 


[C]e  ♦ he([K]e  ♦ pl]e) 


(51) 


All  operations  are  again  limited  to  matrix  inner  Droducts  on  an  elemental 
basis.  The  rank  of  [j]  equals  the  order  of  f<sQ};  specific  (Dirichlet) 
boundary  constraints  are  applied  within  the  evaluation  of  f F } . 
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NUMERICAL  RESULTS 


Test  Case 

As  apparent  from  Table  3,  the  developed  viscous-inviscid  three-dimensional 
interaction  algorithm  for  V/STOL  jet  flowfield  prediction  requires  substantial 
computational  capabilities.  The  immediate  requirement  is  to  assess  the  pivotal 
elements  of  the  algorithm  to  evaluate  its  intrinsic  viability  and  usefulness. 
The  circular  jet,  injected  perpendicular  to  the  plate  and  to  the  crossflow 
represents  the  basic  geometry  and  a substantial  experimental  data  base  exists. 
Therefore,  computational  tests  were  conducted  for  a normal  circular  jet  at 
A = 8 and  dj  = 0.05  m (2  in.).  The  data  of  Mosher  (ref.  10),  also  Fearn  and 
Weston  (ref.  13),  provide  plate  Cp  distributions  for  Step  B comparisons,  and 
the  downstream  interface  shape  for  Step  C can  be  estimated  from  ref.  12,  for 
the  case  of  Uo  = 15  m/s  (50  f/s).  Eiuation  (30)  provides  the  remaining  Step 
A input  requirement  r(xj). 


Ground  Plane  Interface  Contour 

The  first  requirement  of  the  algoritm  is  determination  of  a geometry  for 
the  intersection  of  the  viscous-inviscid  interaction  surface  with  the  injector 
plate.  Step  B of  Table  3.  The  contour  determined  may  not  be  unique,  but  the 
existence  of  at  least  one  contour  is  mandatory.  Viewing  Figure  3,  the  assump- 
tion is  that  the  blockage-only  contour  is  an  elementary  airfoil  belonging  to 
the  Karman-Trefftz  family.  The  adjustment  parameters  are  thickness  ratio  and 
trailing  edge  included  angle,  and  the  airfoil  may  not  close  at  the  trailing 
edge.  The  comparison  plate  pressure  distributions  are  shown  in  Figure  13. 

The  position  of  the  jet  within  the  contour  is  also  an  unknown  to  be  determined. 

The  computational  capability  to  conduct  the  two-dimensional  solutions 
exists  in  the  C0M0C  computer  program  (ref.  22).  Figure  14  illustrates  a 
representative  macro-element  discretization,  with  boundary  condition  specifi- 
cations noted  thereon,  and  the  resultant  finite  element  computational  grid 
generated  by  the  program.  From  Figure  14,  a large  thickness  ratio  appears 
required  for  the  contour,  yet  pressure  recovery  in  the  trailing  edge  region 
must  be  modest.  Figure  15  summarizes  the  results  of  tests,  conducted  with 
large  thickness  ratios  and  large  closed  trailing  edge  angle,  to  serve  as  a 
comparison  basis.  A large  negative  pressure  region  is  induced,  but  trailing 
edge  recovery  is  much  too  large.  The  next  step  assessed  recovery  modification 
by  opening  the  trailing  edge  to  simulate  the  anticipated  blunt  base  requirement. 
Figure  16  summarizes  the  results,  which  indicate  smaller  trailing  edge 
recoveries  can  be  acheived  but  remain  unacceptably  high. 

Reviewing  these  data  indicates  that  a segment  of  the  assumed  airfoil 
contour,  from  the  trailing  edge  forward  to  some  x/C,  must  lie  beyond  the 
interaction  contour  footprint  on  the  plate.  To  ascertain  this,  the  experimental 


Figure  13.  Experimental  Contours  Of  Constant  Surface  Pressure 
Distributions  For  Circular  Jet,  X = 8. 
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Figure  14.  Macro  Element  And  Finite  Element  Discretizations  For 
Karman-Trefftz  Potential  Solutions 
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Figure  15.  Computed  Inviscid  Pressure  Distributions  For  Various  Closed 
Trailing  Edge  Angles 
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pressures  were  interpolated  to  determine  Cp  vs.  x/C  on  various  thickness  ratio 
airfoil  contours,  as  illustrated  in  Figure  13,  assuming  a larger  than  required 


computed  for  various  thickness  ratios  and  trailing  edge  included  angle  and  gaD. 
A summary  of  results  is  shown  in  Figure  17.  The  35%  thick  airfoil,  with  vanish 
ing  trailing  edge  angle  and  gap,  and  with  the  last  35%  of  chord  discarded  met 
the  requirements  of  Step  B.  By  moving  the  center  of  the  jet  just  downstream  of 
its  assumed  initial  coordinate,  Figure  17e,  it  became  possible  to  completely 


envelope  the  computed  with  the  experimental  Cp.  The  downstream  intersection 
of  these  two  curves  determined  the  truncation  point  of  the  airfoil. 


Therefore,  for  this  jet  test  case  at  least,  it  is  possible  to  determine  an 
aerodynamic  blockage-only  contour  that  meets  the  requirement  of  the  interaction 
algorithm.  The  truncated  airfoil  contour  is  completed  by  dropping  a perpen- 
dicular segment  to  the  center  line  axis.  The  local  differences  between  and 
Cp  yield  the  corresponding  entrainment  velocity  distribution  u-1-  using  p 

equation  (4).  Furthermore,  the  jet  has  become  located  within  the  domain. 

Hence,  all  requirements  of  Step  B.,  Table  3 have  been  met. 


Blockage  Plus  Entrainment  Potential  Solution 

Step  C,  the  next  computational  sequence  in  the  algorithm,  has  as  its 
primary  requirement  the  determination  of  the  inviscid  pressure  distribution 
on  the  entire  interaction  surface.  Recalling  equations  (7)-(10),  the  3DPNS 
solution  is  coupled  to  the  3D<}>  solution  via  the  complementary  pressure  distri- 
bution, the  boundary  condition  for  which  is  the  inviscid  pressure.  For  the 
test  case,  the  interaction  surface  contour  at  xx/d  = 0,  after  truncation,  is 
equal  to  3. 2d  and  the  thickness  equals  0.7d. 

The  interface  surface  is  interpolated  wf-th  the  assistance  of  any  available 
data.  Figure  18  illustrates  a three-dimensional  finite  element  discretization 
of  the  3D<f>  solution  domain,  where  the  interface  surface  S^Xj)  has  been 
smoothly  interpolated  over  0 < Xj/d  < 5.2.  The  second  requirement  for  the  3D<t> 
solution  is  distribution  of  the  entrainment  velocity  onto  S(xj).  Experiment 
indicates  it  vanishes  by  Xj/d  = 5.2,  and  in  the  absence  of  other  data,  a linear 
distribution  was  assumed. 

The  computational  requirement  for  Step  C is  solution  of  the  three-dimen- 
sional potential  flow  equation  on  the  domain  shown  in  Figure  18,  subject  to 
the  entrainment  boundary  conditions.  The  basic  capability  to  solve  this  system 
is  resident  within  finite  element  structural  analysis  codes,  and  the  SAP 
program  was  utilized  (ref.  23).  Only  indicative  results  are  obtainable,  how- 
ever, since  the  (displacement)  vector  dependent  variable  character  of  elas- 
ticity requires  "interpretation"  to  run  the  scalar  variable  $.  This  was  accom- 
plished by  equating  <p  to  the  z -component  of  displacement  and  specifying  the 
(x,  y)  plane  entrainment.  The  non-zero  z entrainment  component,  required  for 
Step  E.2  could  not  be  specified. 


1/  V 


a ) t/C  = 0.33,  6 = 20°,  A/C  = 0.16  b)t/C  = 0.38,  S = 20°,  A/C  = 0 


c)  t/C  = 0.25,  6 = 20°,  A/C  =0  d)  t/C  = 0.25,  6 = 40°,  A/C  = 0 


SYMBOL  PRESSURE 


V1- 


Experimental 

Inviscid  Blockage 

Original  Location 


■ / if  W '{-x, 

e)  t/C  = 0.35,  6 = 0°,  A/C  = 0 


Figure  17.  Computed  Pressure  Coefficient  Distributions  On  Various 
Interface  Geometries,  Circular  Jet,  X = 8 


Figure  18.  Three-Dimensional  Potential  Flow  Solution  Domain  For 
Circular  Jet,  A = 8,  UM  * 15m/s 
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Within  these  constraints,  a 3D$  solution  was  obtained  using  SAP  for  the 
specified  interface  contour  and  entrainment  distributions.  Figure  19  illus- 
trates this  solution  with  plotted  as  z-displacements  superimposed  on 
computational  grid  planes.  The  computed  pressure  distribution  on  the  injector 
plane  is  shown  in  Figure  20.  Agreement  with  the  experimental  data  is  quite 
good  in  the  near  vicinity  of  the  interface  contour,  see  Figure  13.  The  influ- 
ence of  three-dimensional  effects  would  be  more  strongly  reflected  away  from 
this  region.  An  accounting  of  the  absent  vertical  entrainment  is  required 
for  completion  of  Step  C.  The  basic  concept  is  verified. 


Turbulent  Jet  Region  Solution 

The  first  requirement  for  the  3DPNS  turbulent  flow  solution  to  be  obtained 
within  the  envelope  of  interaction  surface,  is  initialization  of  all  dependent 
variables.  A finite  element  discretization  of  the  initial  solution  plane  is 
shown  in  Figure  21.  The  previous  step  provides  the  inviscid  pressure  boundary 
conditions  for  solution  of  the  complementary  pressure  distribution.  These 
linear  solutions  are  performed  once  for  each  plane  of  the  3D<j>  solution. 
Quadratic  interpolation  is  then  employed  to  compute  pc  and  the  Xi -derivative 
of  pc  (for  the  u^  momentum  equation)  at  any  xa  station.  Figure  22  is  the 
computed  pc  distribution  at  the  ground  plane,  and  a substantial  transverse 
gradient  has  been  produced.  This  agrees  with  the  observations  of  Fearn  and 
Weston  (ref.  13),  and  will  induce  a substantial  lateral  spreading  perpendicu- 
lar to  the  cross-flow  direction. 

The  location  of  the  jet  is  known  within  the  3DPNS  solution  domain.  The 
initial  distribution  of  uj  inside  the  jet  may  be  specified  non-uniform  depen- 
dent on  available  data.  Since  the  entrained  flow  enters  mostly  behind  the  jet, 
the  non-uniform  distribution  shown  in  Figure  23  was  specified,  with 
Uj  = 1.1  Uj(  = 8.8  Uoo)  in  the  near  region.  Outside  of  the  jet  proper,  the 
axial  velocity  distribution  is  unknown.  However,  elementary  boundary  layer 
concepts  require  that  U!  > 0 to  account  for  displacement  effects.  A uniform 
level  equal  to  0.001  U«,  was  specified. 

As  previously  discussed,  the  inviscid  entrainment  velocity  distributions 
become  boundary  conditions  for  the  3DPNS  solutions  for  transverse  velocities 
uA,  2 <_  i <_  3.  A Uj  distribution  is  required  initialized  throughout  the 
domain  that  matches  the  entrainment  determined  on  Sj  from  Step  B.  One  pro- 
cedure would  be  to  set  the  convection  terms  in  momentum  equations  (15)-(16) 
to  zero  and  solve  the  resultant  two-dimensional  Poisson  equations.  Only  the 
transverse  complementary  pressure  gradients  9pc/ax£  are  known,  however,  since 
the  Reynolds  stresses  involve  the  unknown  turbulence  kinetic  energy  and  dissi- 
pation distributions,  k and  e see  equations  (21).  The  level  of  k is  probably 
small  everywhere  in  comparison  to  that  on  the  periphery  of  the  jet,  however, 
so  a first  order  approximation  is  to  solve  the  Laplacian  for  each  u^  subject 
to  the  entrainment  boundary  conditions.  Figure  24a  shows  the  resultant  dis- 
tribution of  U£  on  Sj  in  the  form  of  a vector  plot.  The  action  of  entrainment 
is  quite  noticeable.  For  comparison,  Figure  24b  is  the  same  solution  with 
iP-  = 0 on  Sj. 
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Figure  22.  Computed  Complementary  Pressure  Distribution  On  3DPNS 
Domain  at  Sj 
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Once  the  mean  flow  velocity  distributions  are  established,  the  initial 
distributions  of  k and  e can  be  determined  using  an  adaptation  of  mixing 
length  theory  to  estimate  local  shear  stress  plus  definition  of  a dissipation 
length  scale.  Recall  that  for  an  elementary  two-dimensional  boundary  layer 
flow  the  dominant  Reynolds  stress  is  evaluated  as 


-U1U2  = 


o2Z2 


(52) 


where  w is  Van  Driest  wall  damping  and  z is  the  mixing  length  (cf. , ref.  20). 
Comparing  equations  (21)  and  (52),  and  generalizing  the  mean  shear  evaluation, 
yields 


C4 


k2 


auj  auj 


(53) 


The  additional  required  relation  to  solve  for  a unique  k and  e is  the 
dissipation  length  scale  z^. 


*d  s 


v e 


3/2 


(54) 


where  Cv  is  a correlation  coefficient.  Using  the  bi-derivative  evaluation 
option  in  C0M0C  III,  assuming  a constant  z ^ that  scales  the  non-dimensional 

computed  k < 1%,  and  assuming  a uniform  background  level  of  ~ .01%  yields  the 
initial  distribution  of  turbulence  kinetic  energy  shown  in  Figure  25.  Notice 
the  non-negligible  levels  essentially  coincide  with  the  jet  periphery. 


This  computational  operation  completes  the  3DPNS  solution  initiation 
sequence  specified  in  Step  D,  Figure  3.  All  required  information  and  data 
are  now  on  hand  for  execution  of  the  3DPNS  solution.  Step  E.  The  entrainment 
aspects  of  the  interaction  algorithm  require  a restatement  of  the  finite 
element  numerical  solution  algorithm  to  account  directly  for  the  continuity 
equation  (13)  in  the  solution  of  equation  (14)  for  Uj.  Equation  (31)  is  the 
expansion  of  the  latter  in  the  cylindrical  coordinate  system.  For  constant 
density,  the  similar  expansion  for  equation  (13)  is 


L(p) 


l_^e  + l 3(rur) 
r 36  r 3r 


(55) 


Both  equations  (31)  and  (55)  contain  the  axial  derivative  of  the  predominant 
velocity  u^  as  an  unknown.  The  transverse  plane  velocity  distributions  u^ 
have  been  initialized  by  the  combined  action  of  blockage  and  entrainment 
with  the  latter  the  mechanism  for  mass  influx  into  the  3DPNS  domain.  The 
associated  added  mass  flux  must  exit  the  3DPNS  domain  only  in  the  direction 
parallel  to  u0;  hence,  equation  (55)  represents  a source  for  ue  and  a 
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Figure  25.  Computed  Initial  Distribution  Of  Turbulence  Kinetic  Energy 


On  3DPNS  Domain 
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differential  constraint  on  the  solution  of  equation  (31).  This  is  accounted 
for  within  the  finite  element  solution  algorithm  in  the  manner  equivalent 
to  use  of  Lagrange  multiplier  a.  Hence,  the  algorithm  statement  of  equation 
(42)  for  u0  becomes 


(N(x^) } L (ue)dt  - a 


{N(x^) } L(p)dx 


-X|{N(xt))  t (ue)  (It 

3R£nsR2 

e 


= {0} 


(56) 


The  first  step  in  the  3DPNS  solution  is  evaluation  of  the  xj  derivatives 
of  the  dependent  variable  vector  {q > . The  Uj  velocity  field  is  of  greatest 
interest,  and  Figures  26-27  are  plots  of  the  derivatives  computed  by  equation 
(56),  for  a = 1 and  a = 0 (no  entrainment).  The  sole  action  in  the  a = 0 case 
is  to  retard  the  jet  core  and  accelerate  the  periphery.  However,  for  a t 0, 
an  acceleration  is  induced  everywhere  by  the  action  of  continuity  in  equation 
(56).  Figure  28  shows  the  Ui  distribution  computed  within  the  jet  region 
after  proceeding  one-half  diameter  downstream,  i.e.  parallel  to  the  xx  axis. 

The  entrainment  and  complementary  pressure  boundary  conditions  previously 
established  in  Steps  B and  C were  employed.  This  action  has  accelerated  the 
Cij  velocities  lying  on  the  interfacial  surface  by  ~ 10%  U*,.  The  required  step 
at  this  station  would  be  comparison  to  the  entrained  30$  solution.  If  agree- 
ment was  satisfactory  the  3DPNS  solution  is  continued.  If  not,  the  algorithm 
would  be  cycled  to  Step  C and  the  3D<(>  solution  recalculated  with  an  adjusted 
exit  velocity  parallel  to  Xj.  This  requires  formal  establishment  of  the 
specific  3D<t>  solution  capability.  For  the  uj  velocity  distribution  in  Figure 
28,  the  increase  in  mass  flux  over  that  at  the  initial  solution  plane  is  11%, 
as  determined  by  integration  of  the  computed  Qi  distributions  over  the  trans- 
verse plane. 
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Figure  26.  Computed  Axial  xx  Derivative  Distribution  Of  Jet  Region  Velocity 


Without  Entrainment 


SUMMARY  & CONCLUSIONS 


A viscous-inviscid  interaction  algorithm  has  been  derived  for  numerical 
determination  of  the  detailed  three-dimensional  flowfield  surrounding  and 
entrained  into  the  wake  of  a V/STOL  jet  issued  orthogonally  into  a subsonic 
crossflow  bounded  on  one  plane  by  an  aerodynamic  surface.  The  dominant 
influence  parameter  coupling  the  three-dimensional  external  potential  flow 
analysis,  and  the  imbedded  three-dimensional  viscous  and  turbulent  analysis 
of  the  jet  and  nearfield  wake,  is  the  entrainment  flow  induced  into  the  jet 
region  and  its  influence  on  the  potential  flow.  The  two  distinct  solution 
domains  are  pressure  coupled  and  the  entrainment-induced  mass  flux  addition 
l to  the  wake  is  the  primary  iteration  variable  of  the  global  algorithm. 


A key  element  of  the  algorithm  is  definition  of  a surface  defining  the 
coincident  boundary  of  the  dual  analyses.  Numerical  studies  have  confirmed 
existence  of  the  interface  surface  for  the  test  case  of  a circular  jet  issued 
into  a crossflow  with  velocity  ratio  equal  to  eight.  This  determination 
yields  definition  of  the  specific  geometry  of  the  interface  intersection  with 
the  injection  plane  and  the  entrainment  velocity  distribution  thereon.  Any 
experimental  surface  pressure  data  may  be  specifically  employed  in  this 
analysis  phase,  and  may  in  addition  be  used  to  complete  definition  of  the 
interface  surface  evolution.  The  summary  three-dimensional  jet  region  compu- 
tations yield  determination  of  the  detailed  evolution  of  the  turbulent  flow 
and  employ  a turbulence  closure  model  involving  all  six  components  of  the 
Reynolds  stress  tensor.  Finite  element  numerics  have  been  employed  to  compu- 
tationally evaluate  the  five  distinct  computational  phases  of  the  algorithm. 


Several  additional  numerical  developments  are  required  completed  to  fully 
assess  the  generality  of  the  developed  interaction  algorithm. 

Specifically: 


1.  A three-dimensional  potential  flow  analysis  with  entrainment  is 
required  established  to  complete  assessment  of  convergence  of  the  dual  analyses 
to  a unique  velocity  field.  The  specific  requirement  is  to  be  able  to  specify 
onset  velocity  distributions  on  all  surfaces  of  the  domain  boundary  and  to 
generate  the  resultant  three-dimensional  potential  flow  distribution. 

2.  Having  acheived  completion  of  the  development  specified  in  1.  above, 
additional  numerical  solutions  for  the  perpendicular  circular  jet  are  required 
to  assess  factors  related  to  convergence. 


3.  The  versatility  of  the  interaction  algorithm  is  required  assessed  by 
application  to  non-perpendicular  injected  circular  jets  and  to  rectangular  jets 
of  practical  aspect  ratio.  The  required  experimental  data  on  (primarily) 
surface  pressure  distributions  are  available. 

4.  The  compromising  influence  induced  into  the  three-dimensional  turbu- 
lent jet  region  analysis  by  the  imposed  parabolic  approximation  should  be 
assessed.  This  is  primarily  manifested  in  the  computed  entrainment  velocity 
distribution  on  the  injection  plane  and  in  the  resultant  nearfield  induced 
pressure  distribution. 

5.  The  developed  interaction  algorithm  should  be  critically  evaluated  for 
potential  simplification  to  yield  an  entrainment  model  applicable  to  a compu- 
tational ly-elementary  strictly  potential  flow  analysis. 
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APPENDIX 


THE  THREE-DIMENSIONAL  PARABOLIC 
NAVIER-STOKES  EQUATIONS 


The  description  of  a state  point  in  multi-dimensional  fluid  mechanics  is 
contained  within  solution  of  a coupled  nonlinear  partial  differential  equation 
system  describing  conservation  of  mass,  momentum  and  energy.  Unique  solutions 
are  obtained  upon  closure  by  specification  of  constitutive  behavior  and 
boundary  conditions.  In  Cartesian  tensor  notation,  the  non-dimensional  con- 
servation form  of  the  Navier-Stokes  equation  system  is 


Up) 

Ucui>  = * siT 

JpujUj  + p6.jj  - Re“1a1  + ip^-  = o 

L(pH)  3 ^<pH)  * ,£- 

3 

pujH  - (y- l)M2Re_1aijUi  - Re^Pr^u^- 

J, 

- (Y-weH-o 


(A.  1) 

(A. 2) 

(A. 3) 


The  dependent  variable  in  equations  (A.1)-(A.3)  have  their  usual  inter- 
pretation in  fluid  mechanics,  i.e.,  p is  mass  density,  Uj  is  the  velocity 
vector,  p is  the  static  pressure,  b is  a body  force,  and  H is  the  stagnation 
enthalpy.  Furthermore,  o-jj  is  the  Stokes  stress  tensor,  defined  as 


a.  • = u 
U M 


+ % 

3xj  8xi 


2w  3uk 
T3x“6ij 


(A. 4) 


where  u is  the  dynamic  viscosity.  The  equation  of  state  for  a perfect  gas 
closes  the  definition. 

Equations  (A.1)-(A.4)  are  valid  for  both  laminar  and  turbulent  flows. 
Their  solution  for  the  latter  becomes  tractible  in  a practical  sense  only 
after  time-averaging.  Mass-weighted  time-averaging  is  employed  and  the 
Reynolds  decomposition  is  (ref.  A.l) 

Ujlx,,*)  = a J (X, ) ♦ Uj'(x(,t)  ( 

The  mass-weighted  time-average  velocity  is  defined  as 


Sj  5 pUj/5 


(A.6) 
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and 


puT  = lim 


t+T 


(pui  - pu^)dt  = 0 


This  operation  yields  the  important  relations 


pu.u  . = pu.u . + pu.u . 
i J i J i J 

H = h + J a.i.  + £ puruT/p 


(A. 7) 

(A. 8) 
(A. 9) 


Substitution  of  equations  (A.5)-(A.6)  into  (A.1)-(A.4),  time-averaging 
and  collecting  terms  yields  the  time-averaged  Navier-Stokes  equation  system 
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(A. 12) 


The  three-dimensional  parabolic  Navier-Stokes  (3DPNS0  equations  are  established 
from  equations  (A. 10)- (A. 12)  under  the  assumptions  that: 


(1)  A predominant  mean  flow  direction  is  uniformly  discernible. 

(2)  In  this  direction  (only)  diffusion  processes  are  negligible 
compared  to  convection,  and, 

(3)  Overall  three-dimensional  elliptic  character  is  provided  by 
interaction  with  the  potential  freestream  flow. 


The  parabolic  approximation  to  equations  (A. 1 1 ) - ( A . 12)  basically  con- 
strains summation  limits.  Assuming  the  xj  coordinate  aligned  with  the 
direction  of  predominant  flow,  the  parabolic  approximation  is  concisely 
expressed  as 


'5.  ♦ . 1 * 

3Xj  3x.j  3 3xk  °ij 


(A. 13) 
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The  subscript  bar  notation  denotes  the  index  not  eligible  for  the  summation 
convention,  but  identical  with  the  value  of  the  synonymous  tensor  index. 

The  parabolic  approximation  to  equation  (A. 12)  also  requires  eliminating  the 
xj  component  of  the  diffusion  terms  involving  the  time-averaged  viscosity  ii. 
The  parabolization  procedure  effectively  removes  the  boundary  value  character 
on  the  xj  coordinate  from  the  equations. 

The  3DPNS  equation  system  becomes  closed  upon  identification  of  the 
components  of  the  Reynolds  stress  tensor.  Baker  (ref.  A. 2)  has  developed  a 
constitutive  equation  for  the  kinematic  Reynolds  stress  tensor  - u^uj  as 
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ij 
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au.  au.' 
— L + — J. 
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au.  au. 
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’Mk  + 
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(A. 14) 


This  expansion  results  from  reexpression  of  triple  correlations  within  the 
Reynold's  stress  transport  equation  using  the  model  of  Launder,  Reece  and 
Rodi  (ref.  A. 3),  and  is  a generalization  of  the  original  analysis  by  Gessner 
and  Emery  (ref.  A. 4).  In  equation  (A. 14)  ajj  is  the  diagonal  tensor 


aij  * 3k  (u£ujp  ai6ij 


(A. 15) 


and  equation  (A. 14)  represents  an  anisotropic  tensor.  In  equation  (A. 15)  the 
aj  are  the  coefficients  that  admit  anisotropy  and  ax  = a2  = a3  = C3  is 
suggested,  where  the  C-j  are  functions  in  terms  of  two  "universal"  empirical 
constants  derived  by  Hanjalic  and  Launder  (ref.  A. 5).  Additionally,  k is  the 
kinetic  energy  of  the  turbulence  velocity  field 


k 5 

1 

and  e is  the  isotropic  dissipation  function  defined  by  the  contraction 


(A. 16) 
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The  governing  differential  equations  for  k and  e are  (ref.  A.l) 


(A. 17) 
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The  various  C|<  and  Ca  are  empirical  model  constants. 


The  3DPNS  equation  system  is  obtained  by  deletion  of  the  time  derivatives 
in  equations  (A. 10)-(A. 11)  and  (A. 18)-(A. 19)  and  eliminatinq  the  x:  stress 
components  in  equation  (A. 14)  in  the  manner  illustrated  for  the  Stoke's  stress 
tensor.  The  resultant  3DPNS  equation  system  for  isoenergetic  flow  is 
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(A. 20) 
(A. 21) 

(A. 22) 
(A. 23) 
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These  equations  introduce  the  3DPNS  limited  summation  index  convention 
1 <_  i,  j £ 3,  2 £ k,  i <_  3.  The  di  1 i tat  ion  term  in  the  Stoke's  stress 
tensor  has  been  deleted  from  equations  (A.21)-(A.23)  as  neglibibly  small. 

The  parabolic  approximation  to  equation  (A. 14)  is  used  in  completing  terms 
in  equations  (A.21)-(A.25),  yielding 
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For  equation  (A. 21),  the  divergence  term  becomes 
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Defining  the  "effective"  kinematic  viscosity 


v9  4+C4  F k 


equation  (A. 28)  becomes 


(A. 28) 
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In  the  transverse  flow  momentum  equations  (A.22)-(A.23),  the  divergence 
terms  become 
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(A. 30) 
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The  first  right-hand  term  is  basically  identical  to  equation  (A. 28)  for  Uj, 
with  the  addition  of  terms  involving  shear  of  the  alternative  transverse 
plane  mean  velocity  component  u£.  The  second  term  is  a non-homogeneity 
involving  shear  components  of  the  predominant  velocity  Uj  which  acts  as  a 
source  term  within  the  transverse  momentum  equations. 


The  divergence  term  common  to  both  the  k and  e equations  becomes. 


J_  2C  5 *1-23. 

3x£  3UqP  e 3x£ 


(A. 31) 


where  q represents  eigher  k or  e,  and  is  identical  to  the  Reynold's  stress 
contribution  in  equation  (A. 28)  with  Ci,  replaced  bv  C_.  The  production  terms 
for  both  k and  e are  basically  identical,  with  the  latter  multiplied  by  C*  k . 

The  Reynold's  stress-velocity  shear  contraction  under  the  parabolic  approxima- 
tion becomes 


— v k2 

fVk  * pc»  T 


3Qj  aQj 


+ 2pC2  |K  ‘ k 


K1 

3Qj  3Uj 

3xx 

(A. 32) 


The  key  element  of  a viscous-inviscid  interaction  algorithm  is  the  pressure 
coupling  between  the  two  flowfield  descriptions.  Derivation  of  the  pressure 
coupling  algorithm  employs  the  steady  parabolized  mean  flow  momentum  equations 
(A. 21)- (A. 23),  and  is  based  on  the  two-dimensional  transverse  plane  Poisson 
equation  obtained  by  application  of  the  divergence  onerator.  Recalling  the 
limited  summation  convention,  1 < j < 3,  and  2 <_  k,  a <_  3,,  this  equation  is 


L ( p) 


iL  + 

Sx£  3Xj 


(5V*> 


(A. 33) 


The  algorithm  is  based  upon  the  observation  that  the  solution  of  a linear, 
elliptic  Poisson  equation  consists  of  a complementary  and  a particular  con- 
tribution as 
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(A. 34) 


P(xi ) = Pc<xi>  + Pp(xi) 


By  definition,  the  complementary  solution  satisfies  the  homogeneous  form 
of  equation  (A. 33),  i.e. 

32pc 

L(pc)='^|~=0  (A. 35) 

subject  to  the  known  boundary  conditions  for  p(x.,-).  Since  the  bounding 
inviscid  flow  pressure  distribution  is  assumed  everywhere  known, 

P(*t)  = pI(*4)  (A. 36) 


In  equation  (A. 36),  x^  indicates  xz  constrained  to  the  boundary  of  the 

viscous  solution  domain,  9ft  = 9R2  * xj,  and  pi  is  the  (inviscid  flow)  pressure 
level  distribution.  Elsewhere  on  9ft,  the  appropriate  boundary  condition  is 
vanishing  normal  gradient,  i.e. 

£(pc)  = 9^  pc  * ni  = 0 (A. 37) 

The  particular  solution  pp  is  any  function  satisfying  equation  (A. 33) 

with  homogeneous  boundary  conditions  on  closure  segments  of  ft  coincident  with 
p1  known.  Elsewhere,  the  boundary  condition  for  pp  is  established  from  the 

inner  product  of  the  parent  steady  momentum  equations  (A. 22)-(A. 23)  with  the 
closure  surface  normal.  Since  the  convection  term  vanishes  identically  at  a 
wall  or  symmetry  plane. 


£(pp}  + ^ "l  ' 3^  [>  ‘ pukX]n*  = ° (A-38) 

The  appropriate  Poisson  equation  for  pp  is  the  inviscid  form,  (ref.  A. 6) 


L(pp) 


_9_ 

3x, 


liu  jl 

8xt  3xj 


““A  * puA-j 


= 0 


(A. 39) 


Equations  (A. 33)- (A. 39)  define  the  algorithm  for  determination  of  pressure 
distributions  in  the  three-dimensional  parabolic  viscous  flow  solution.  The 
*i  pressure  gradient  used  in  solution  of  the  predominant  flow  momentum  equation 
is  determined  solely  from  the  complementary  pressure  solution.  The  complemen- 
tary and  particular  solutions  are  summed,  as  in  the  defining  equation,  and 
differentiated  to  form  the  transverse  pressure  gradient  field.  The  continuity 
equation  (A. 20)  is  substituted  directly  into  the  right  side  of  equation  (A. 39). 
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